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© A process and an apparatus for measuring the biocurrent distribution are disclosed. Multiple current more 
than the number of isolated currents possibly existing in a living body are assumed and the virtual magnetic field 
created by the assumed currents at actual measuring points (12-1, .... 12-M) is calculated. The estimated 
location and current vectors of said assumed currents are consecutively changed to decrease the difference 
between said virtual and actual magnetic fields so that the estimated location and current vectors of the 
assumed currents having the number corresponding to the true isolated currents substantially approximate to the 
true values and the estimates of the other assumed currents approximate to zero. This enables the biocurrent 
distribution to be exactly determined even when the number of isolated currents (q n ) is not preliminarily known. 
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BACKGROUND OF THE INVENTION 

The present invention relates to a process and an apparatus for measuring the biocurrent distribution 
and in particular to a process and an apparatus for indirectly measuring the biocurrent distribution by 
5 measuring the biomagnetic field and estimating the distribution of the biocurrents which are sources of the 
magnetic field based upon the measured data. 

Estimation of the biocurrent distribution from the measured values of the biomagnetic field is usually 
performed as follows: 

Description will be made by way of a neuromagnetic measurement which is a typical biomagnetic field 
10 measurement. 

The coordinate system is defined as shown in Fig. 1. It is assumed that a current represented as a 
vector q n exist in position r n in the drawing. Such an isolated current vector can be assumed as far as the 
interface of the brain can be considered as spherical surface in the field of neuromagnetic field measure- 
ment and this current vector will be referred to as current dipole. A reference r m in Fig. 1 denotes the 
15 coordinate of measuring point. 

If it is assumed that N current dipoles exist in the brain, the magnetic flux density vector B m at point r m 
can be expressed from Biot-Savart's law as follows: 



20 N 

S 

n=l 4 * I r m -r 
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„ „ MO qn * (r m - r n ) . 
Bm = £ ; 1 1 ) 



25 The magnetic field component which can be measured by a magnetometer is the component of B m 
normalto the surface of the brain. The actual value measured at r m is denoted as D m . 

If estimated location and current vectors of each dipole is expressed by superscripting a reference mark 
A , the virtual measured value D m which is calculated from the estimates is expressed as follows: 



N 

n=l 4 * I r m - t n |3 I ml 



- n y ° qn* (r m -r n ) « r m 
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A cost function is herein defined as follows: 

r 2 r ... ? N ; qi/ q2/ <3n) 



M 

= 2 (Dm - D m )2 (3) 

m=l 



Equation (3) shows the degree of matching between an estimated data and an actually measured data. 

Optimum estimated values have heretofore been determined as ?i . .... ? N , qi q N which minimize the 

cost function defined by the equation (3). 
so It is found from the foregoing that the number N of the biomagnetic field sources, that is, the current 
dipoles should be known in order to carry out conventional measuring method. In other words, the 
conventional method can be carried out only if the number of the current dipoles is known. 

Jf N which is different from the number of actual current dipoles is preset and a solution of n , .... ? N , q, , 
.... qw which minimizes the cost function is determined, the solution would then show a distribution different 
55 from the actual biocurrent distribution. 

SUMMARY OF THE INVENTION 
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It is an object of the present invention to provide a process and an apparatus for measuring the 
biocurrent distribution which is capable of obtaining a solution exactly showing the biocurrent distribution 
even if the number of the magnetic field sources, that is, the exact number of the current dipoles in a living 
body to be measured is not known. 

5 In order to accomplish the above-mentioned object, in accordance with the present invention, multiple 
currents larger than the numbers of the isolated currents which are expected to exist in a living body are 
assumed, a virtual magnetic field created by the assumed current at the actually measured point is 
calculated. Estimated location and current vectors of the assumed current are changed to reduce the 
difference between the virtual magnetic field and the actual magnetic field for determining the optimum 

w estimates. In this case, the absolute values of the excess assumed current vectors are approximated to 
zero. 

If the multiple isolated currents, for example, the current dipoles larger than the number of currents 
actually existing are assumed and the estimates are changed so that the matching between the distribution 
of the magnetic field calculated from the assumed current dipoles and the distribution of the actually 

75 measured magnetic field is enhanced, an excessive freedom degree would be introduced. Accordingly, 
there is the high possibility that a solution different from the actual current distribution is obtained. In other 
words, it is to suppose that the magnetic field formed at the measuring point may be substantially equal to 
the magnetic field formed by the current dipole assumed as true current dipole (refer to equation (1)), if the 
current vector of the current dipole which is assumed to locate in a position closer to the measuring point 

20 than the position of the true current dipole is smaller than the true current dipole. 

Hence, in accordance with the present invention, a fact is considered that the assumed current dipoles 
are more than the true current dipoles in magnitude by the number of the assumed current dipoles over the 
number of true current vector of the assumed current dipoles if the assumed current dipoles are compared 
with the true current dipoles in only magnitude or the energy (for example, magnetic field) which they give. 

25 That is, when optimum values of the estimated locations and current vectors of the currents are 
determined so that the difference between the magnetic field calculated from the assumed currents and the 
actual magnetic field becomes smaller, the excessive freedom degree of the optimum values is restricted 
by decreasing the magnitude of current dipoles or the energy that they give so that the optimum values of 
the current vectors of the assumed currents which is an excess over the number of the true currents 

30 approximate to zero. 

In the present invention, an excessive freedom degree is restricted by using the total sum or the total 
sum of a-th power of the absolute values of the current vectors of the assumed current dipoles, or the total 
sum of the magnitudes of the magnetic fields, each formed solely by each assumed current dipole or the 
total sum of a-th power of the magnitudes thereof. 

35 

BRIEF DESCRIPTION OF THE DRAWINGS 

Fig. 1 is a schematic view illustrating an embodiment of neuromagnetic field measuring of the present 
invention and the coordinate system used for estimation of magnetic field sources; and 
40 Figs. 2A and 2B are flow charts showing a process of simulated annealing of the present embodiment. 

DESCRIPTION OF THE PREFERRED EMBODIMENTS 

Embodiments of the present invention will be described with reference to a neuromagnetic field 
45 measurement. 

The neuromagnetic field measurement is conducted by measuring the component of the magnetic field 
in a direction normal to the brain surface at measuring points around a head while moving a magnetic flux 
meter 10 from points 12-1 to 12-M as shown in Fig. 1. Presumption of space coordinates and current 
vectors of current dipoles in the brain which are sources of generating magnetic fields can be conducted 
so based upon the measured values at measuring points around the head is accordance with present invention 
as follows: 

The number of the measuring points on the surface of the head is assumed as M and the measuring 
points are numbered from I to M. The position of the m-th measuring point is represented by a vector r m . 
The measured value of the magnetic field at that point Is represented as D m . 
65 It is assumed that N current dipoles are contained in the brain. The present invention is applicable to a 
case in which an exact number of the dipoles are not know. The N is for a time preset larger than the 
number of the dipoles which are expected to exist in the brain. 

The estimated location and current vectors of the thus assumed current dipole are represented ? n and 
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q n , respectively wherein n = 1, 2, .... N. The magnetic field D m (m = 1, 2 M) which is created at the 

measuring points by the thus assumed current dipoles is determined by an calculation. That is, D m can be 
referred to as the measured value of a virtual magnetic field determined from the assumed current dipoles. 
If the surface of the head is approximated by a spherical face and the center of the sphere is assumed as 
an origin of the coordinates, 6 m is calculated by an equation as follows: 



n=l 4 * I r m -r n l 3 I r m I 



By an optimization operation which will be described in detail, the estimated location and current 
rs vectors of the assumed current dipoles are changed so that a minimum value which the cost function 
assumes is determined. 

Following equation is used as the cost function in the present embodiment. 

20 E(r 17 r 2 , r N ; q ± , q 2 , . # ., § N ) 

M 

25 * £ < D m ~ D m) 2 + Es(q ir q 2 , q N ) 

m=l 



(4) 



30 

The first term of the right formula of equation (4) denotes matching between the distribution of the 
actually measured magnetic field and the distribution of the magnetic field created by the assumed current 
dipoles. If the respective distributions are represented by row vectors G and G, the first term can be 
represented as |G - G| 2 wherein 
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Dm 




On the other hand, the second term in equation (4) is assumed as the sum total of the magnetic fields, 
each created by individual assumed current dipole. The detailed form of Es (qi , q 2 .... q N ) is as follows: 



N 

Es(q lf q 2 , . . . f q N ) = w E I g n I (5) 

n=l 



wherein 
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9n = 



/ 



V 



and g m n is calcutated as follows: 



9 2 n 



9 M r 



(6) 



/5 



9 m n 



Po q * ( r m - r n ) » r m 

4tt lr m -f n l 3 • lr m l 



(7) 



20 



The component g m n in a direction normal to the surface of the brain (the component vertical to the face of a 
detecting coil) of the magnetic field created at the m-th measuring point by n-th current dipole is calculated 
by equation (7) for each m = 1 , 2 M. The magnitude of the row vector g n having respective components, 

25 that is, the norm |g n | of the magnetic field created by the n-th current dipole is calculated. Specifically, the 
root of the sum of the square of the normal component of the magnetic field created by n-th current dipole 
at each measuring point is calculated. The magnetic field |q n | individually created by each assumed current 
dipole is determined in such a manner and the sum total with respect to n assumed current dipoles is 
determined to provide Es in equation (5). 

30 Es (qi q N ) represented by equation (5) serves to suppresses the excessive freedom degree 

introduced by assuming the dipoles move than the number of actually contained dipoles. Although it is hard 
to logically and exactly verify this in case of non-linear inverse problem like the estimation of the magnetic 
field source in the brain which in an object of the present invention, linear inverse problem can be 
mathematically verified. This verification is discussed In H. Haneishi et al. "An analysis of cost function used 

35 in simulated annealing for CT image reconstruction". Applied Optics in Press. 

"w B in equation (5) is a constant for weighting the first and second terms in the right expression of the 
equation (4). ,? w n is preset by preliminarily calculating AE as a trial so that the contribution Is not biased to 
only one of the first and second terms in equation (4) when a change in the cost function AE is calculated 
for optimization operation which will be described hereafter. 

40 Although equation (5) may be directly used as term of Es in equation (4), it is more preferable to use as 
reformed equation as follows: 



N 

45 Es(q ir q 2 , q N )= w 2 | q n | a (5)' 

n=l 

wherein an index of power a is a positive real number. 
The role of a in equation (5') is as follows: 
so True current distribution is estimated by assuming the current dipoles more than the number of the 
current dipoles which are expected to exist in the brain in the present Invention. Accordingly, there is the 
possibility that a plurality of current dipoles are estimated in substantially same position, a in the equation 
(5') determines the nature of the solution obtained in such a case. A case in which two current dipoles u-th 
and v-th are estimated in substantially same position, that is, a case in which ? u = ? v will now be 
55 exemplarily described. 

(1) In case of 0 < o < 1, either one of the current vectors q u and q v of the u-th and v-th current dipoles 
becomes very small. That is, lq u l == 0 or |q v | = 0. In other words, one current vector of a plurality of 
current vectors having interposed space coordinates approximates to the current vector of a current 
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dipole which inherently exists in that position and the magnitudes of the current vectors of the other 
current dipoles become substantially zero. 

(2) In case of a > 1 , q u == q v - These values are about a half of the value of the current vector of the 
current dipole which inherently exists in position r u = r v , 
5 (3) In case of « = 1 , a value q u + q v approximates to a current vector value of the current dipole which 
exists in position r u = r v . Which value a ratio lq u | to lq v l assumes cannot be predicted. 
Briefly, it is preferable that a be other than 1, and 0 < a < 1 is most preferable (for example, a = 1/2). 
An optimization operation which determines estimates ? n . Qn (n = 1, N) for making minimum the cost 
function E (?i. .... ? N ; cji, q N ) defined by equation (4), more specifically equation (5) or (5*) will be 
10 described in detail. The cost function E is non-linear with respect to the variable and has a localized 
minimum value. Accordingly, it is impossible to determine a true minimum value (a global minimum value) 
by known non-linear optimization approaches. The present embodiment adopts a recently proposed 
simulated annealing method which can determine a global minimum value even in such a case. 
This method will now be described with reference to Figs. 2 A and B. 
75 The above-mentioned N current dipoles are assumed and a proper value is preset as an initial value for 

respective coordinate r n (n = 1, .... N) and current vector q n (n = 1 N) is preset as a variable X n (n = 1, 2 

.... N) and the current vector q n (n = 1, 2 N) of each current dipole is preset as a variable X n - 

(n = N + l,N + 2, .... 2N). The values of the cost function E which has been described in detail are calculated 
by using the actually measured and recorded values D m (m = 1, M) of the magnetic field (step 101). 
20 Thereafter, the value of X n is changed by a very small amount (this is referred to as trial). Whether or 
not the trial is accepted is determined with reference to a change AE in the cost function E occurred by the 
trial. The value of each X n is gradually approximated to the true value. 

Firstly, a sufficiently high parameter T is preset. The parameter T determines the probability that the 
change is accepted even if AE > 0, and is referred to as temperature. 
25 AE is calculated by a trial calculation and T isjpreset ^so that exp (-AE/T) = 0.8 to 0.9. The changes |Ar| 

and |Aq| in one trial are preset for variables ?i r N and qi, q 2 , 5n. respectively. (stepJ02) 

It is adequate that the value of |A?| be about 1 mm to 10 mm and the value of |Aq| be about 1/10 to 
1/1000 of a predicted maximum value of |q|. 

Then, respective initial values zero are preset for parameters N T , Ni, N 2 for repeated control of the trial 
30 (step 103). 

One variable X n which will be an object of the trial is then chosen (step 106) to determine a change in 
vector AX n (step 107). If the variable X n is a coordinate, following equations are used. 



35 



40 



(Ax n )x = I Ar Isin 9 cos t|> 



(Ax n )y = I Af Isin 9 cos ip V ( 8 ) 
(A x n )z = I A f I cos 9 



J 

45 If the variable X n is a current vector, the following equations are used, 



50 



55 



(Ax n ) x = I Aqlsin 9 cos 

(Ax n )y = I Aqlsin 9 cos V ( 9 ) 

(Ax n ) z = I Aq I cos 9 

The change in vector is determined by generating uniform random numbers 6 and ^ in the ranges of 0 
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£ $ < 2tt and 0 £ </> £ -n-, respectively. 

The variable X n is replaced with Xn + AX n in step 108. In other words, a trial is executed. 

In step 109, the difference AE between previous and subsequent cost functions E which give a change 
AX n is calculated as follows: 

5 

AE = E(xi, x n + Ax n x 2N ) -E(xi x n x 2N ) (10) 

If AE < 0, the change AX n is then accepted in step 111. That is, calculation program is returned to step 106 
for trial of next variable X n+1 via steps 112 and 105 while X n + AX n is taken as X„. If AE > 0, whether the 

w change AX n is accepted or refused is determined according to a probability P(AE) = exp (-AE/T) 
depending upon temperature as shown in steps 113, 114 and 116. If it is determined that the change is 
accepted, the calculation program returns to step 106 via steps 115 and 105 while taking X n + AX„ as new 
X n . If the change is determined to be refused, X n and the cost function E are returned to original values In 
step 117 and the program is returned to step 106 via steps 118 and 105. 

rs Such a series of steps are executed for each variable, that is, for n = 1 through 2N. This operation is 
repeated total N T ma * times. It is adequate that, for example, a value of about 100 to 400 is adequate as 
N T max . 

Ni in step 1 12 of a flow chart of Fig. 2 A denotes a number of changes accepted in such a direction that 
the cost function is decreased, N2 in step 115 denotes the number of changes accepted in such a direction 

20 that the cost function is increased. N3 in step 105 denotes the number of the refused changes. When it is 
found in step 105 that N T max trials are completed, the program proceeds to step at which determination 
whether or not the temperature T determining the probability of acceptance of change is lowered is 
performed. That is, it is determined based on a value of |Ni - N 2 | whether or not the number N1 of 
increases of the cost function of the accepted trials approximates to the number of decreases of the cost 

25 function. If this value is not less than e (for example, about 0.02), the program is returned to step 103 
without changing the temperature and trials will be repeated N T max times. If the value of IN1 - N 2 |/N is lower 
than e, the program proceeds to step 120. By lowering the temperature T, the probability that the change to 
make AE > 0 is accepted is decreased and thereafter the program proceeds to step 103 and the trials are 
similarly repeated N T max times. 

30 An approach to make T k = $T k _ 1t as k-th temperature step is adopted in the present embodiment to 
lower the temperature T in step 120. A value of about 0.9 through 0.95 is used as * at this time. Another 
approach to lower the temperature by T k = TJ(1 +k) or T k = Vlog (e + k) is also proposed in a literature A 
Vol. 122, p157, 1987 (H. Szu et al., "Fast Simulated Annealing", Physics Letters A, Vol. 122, p157, 1987. 
As N T max trials of changes of each variable and repeated in some cycles while gradually lowering the 

35 temperature T as mentioned above, each variable approximates to an optimum value and the accepted 
change is gradually decreased. When it is determined that N1 + N 2 = 0, in step 118, that is, there is no 
accepted change of the tried changes, a variable which minimizes the cost function, that is, optimum 
location and current vectors of a current dipole are obtained. Therefore, simulated annealing operation is 
ended. 

40 Although the present embodiment has been described with reference to a case in which simulated 
annealing is used for optimization operation, the present invention is not limited to this case. The other 
algorithms which can find out a true minimum value from a function containing a localized minimum value, 
such as Genetic Algorithm proposed by J. H. Holland can be used. 

If areas at which di poles exist are a priory known in advance, the priory information can be incorporated 
45 into the optimization operation as follows: 

In this case, each estimate which minimizes E which is calculated by the following equation (11) in lieu 
of equation (4) is determined. 
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E(r^ f •••/ r N ; q^ f 0$) 



M 

= 2 (Dm - D m )2 + Esfq^ q N ) 
m=l 

70 

+ E L (r lf r N ) (11) 



75 E L is defined as follows: 

1) If all?i ? N exist in preliminarily preset areas, E L (?i, .... r N ) = 0 ^ 

2) If any of ?i, .... ? N moves out of the preliminarily preset areas, E L (ri, .... ?n) - A wherein is a value 
which is sufficiently high for the current temperature Ti that is, such a value that 



20 



25 



exp <-A/T)*0. 

exp By defining E L in such a manner, the probability that a change by which each ? n moves out of the 
preliminarily preset area is accepted Is approximately zero. An estimate which minimizes 



30 M 

= E (D* - D»)2 + Es(§ ir ^) 
m*l 

35 can be determined in preliminarily preset area. 

Although a case has been described in which the second amount expressed by equation (5) or (5'), that 
is, the magnetic field created by the assumed current dipole is used as the second term of equation (4) for 
restricting the freedom degree, equation (5") of the total sum of absolute values of current vectors of the 
presumed current dipoles or equation (5'") of the a-th power of the total sum may be used in lieu of 

40 equation (5) or (5'). 

N 

Es(§ lf q 2 , q*) = w E I qn I (5") 

45 n=l 

N 

Es(q x , q 2 , §m) » w 2 I qn 1° (5'") 

50 n«l 

In equation (5") and (5"*), the meaning of w and a are same as those of equations (5) and (5'). 

Since terms of position vectors are included in equations (5) and (5"), the freedom degree is more 
55 restricted than in equations (5 f *) and (5"'). 

Claims 
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1. A process for measuring the distribution of the biocurrent by measuring the biomagnetic field generated 
by the activity of a living body at a plurality of measuring points (12-1, .... 12-M) to indirectly determine 
the distribution of the isolated currents (q n ) in said living body from the plural measured values of the 
magnetic field, comprising the steps of: 

5 measuring the actual magnetic field at said measuring points (12-1 12-M); 

assuming multiple currents move than the number of said isolated currents possibly existing in said 
living body to preset the locations and the current vectors of the assumed currents; 

calculating a virtual magnetic field at each measuring point from the estimated location and current 
vectors; and 

10 determining including consecutively changing the estimated location and current vectors of said 

assumed currents to decrease the difference between said virtual and actual magnetic fields so that the 
estimated location and current vectors of the assumed currents having the number corresponding to 
the true isolated currents substantially approximate to the true values and the estimates of the others 
assumed currents approximate to zero. 

75 

2. A process for measuring the biocurrent distribution as defined in Claim 1 in which said determining 
step includes the steps of: 

a changing said estimated location and current vectors; 

calculating a first amount representative of the difference between said virtual and actual magnetic 
20 fields; 

calculating a second amount from each of said estimated current vectors; and 

determining said estimated location and current vectors from said first and second amounts. 

3. A process for measuring the biocurrent distribution as defined in Claim 2 in which said step of 
25 determining said estimated location and current vectors includes presetting a cost function including at 

least said first and second amounts and determining the estimated location and current vectors by 
using the cost function. 

4. A process for measuring biocurrent distribution as defined in Claim 3 in which said cost function 
30 includes a linear sum of said first and second amounts. 

5. A process for measuring the biocurrent distribution as defined in Claim 4 In which said step of 
determining said estimated location and current vectors is executed by minimizing the cost function. 

35 6. A process for measuring the biocurrent distribution as defined in Claim 4 in which said second amount 
is the total sum of the absolute values of said estimated current vectors. 

7. A process for measuring the biocurrent distribution as defined in Claim 4 in which said second amount 
is the total sum of the a-th power (a is a positive real number) of the absolute values of said estimated 

40 current vectors. 

8. A process for measuring the biocurrent distribution as defined in Claim 4 in which said second amount 
is the total sum of the absolute values of the magnetic fields, each created at each measuring point by 
each of said assumed current vectors. 

45 

9. A process for measuring the biocurrent distribution as defined in Claim 4 in which said second amount 
is the total sum of the absolute values of the magnetic field's components normal to the surface of said 
living body, each magnetic field being created by each of said assumed current vector at each 
measuring point. 

50 

10. A process for measuring the biocurrent distribution as defined in Claim 4 in which said second amount 
is the total sum of a-th power (a is a positive real number) of absolute values of the magnetic field's 
components normal to the surface of said living body, each magnetic field being created by each of 
said assumed current vector at each measuring point. 

55 

11. A process for measuring the biocurrent distribution as defined in Claim 7 or 12, in which a + 1, 
preferably 0 < a < 1 . 
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12. A process for measuring the biocurrent distribution as defined in Claim 5 in which said step of 
determining said estimated location and the current vectors Includes a first step of accepting a change 
in said estimates executed in the step of changing said estimates when said cost function decreases, 
and a second step of accepting a change in said estimates in a probability relating to a given 

5 temperature even when said cost function increases. 

13. An apparatus of measuring the biocurrent distribution by measuring the magnetic field generated by the 
activity of a living body at a plurality of measuring points (12-1, .... 12-M) to indirectly determine the 
distribution of isolated currents (q n ) in the living body from the plural measured values of the magnetic 

io field, comprising: 

means for measuring the actual magnetic field at said measuring points (12-1 .... 12-M); 
means for assuming multiple currents more than the number of said isolated currents possibly 
existing in said living body to preset the location and the current vectors of the assumed current; 

means for calculating a virtual magnetic field at each measuring point from the estimated location 
15 and current vectors; and 

determining means including consecutively changing the estimated location and current vectors of 
said assumed currents to decrease the difference between said virtual and actual magnetic fields so 
that the estimated location and current vectors of the assumed currents having the number correspond- 
ing to the true isolated currents substantially approximate to the true values and the estimates of the 
20 other assumed currents approximate to zero. 
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FIG. I 




® ® 



FIG. 2B 

a" 




CALCULATE 
P=exp(-AE/T) 



GENERATE UNIFORM 
RANDOM NUMBER R 
OF 0< R<l 



YES 




Xn~*~ Xn~AXn 
E — E ; 



N3—N3+I 
Nt*Nt+I 



13 



14 



-II7 



18 



11 



EP 0 449 231 A1 



FIG. 2A 
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